/**
 * visSteamlineClustering.cpp 
 * Linux(ubuntu8.04).
 * Copyright Thss, THU
 * Author : luxinlai
 */

#include "visStreamlineClustering.h"
#include <ctime>

namespace visualization
{
/****************************** COMPUTE STREAMLINES *******************************/
void visStreamlineClustering::computeStreamlines(void)
{	
//  /******************* Set sources *******************************/
//  // Randomly generate m_sourceNumber numbers for x, y and z respectively.
//  Point *sources = new Point[m_sourceNumber];
//  double x, y, z;
////	cout << "Currently the streamlines are not generated randomly..." << endl;
//  srand( (unsigned int)time(NULL) );
//  for(unsigned int i = 0; i < m_sourceNumber; i++)
//  {
//    x = rand()/double(RAND_MAX);
//    y = rand()/double(RAND_MAX);
//    z = rand()/double(RAND_MAX);
//    
//    x = m_dataField -> getXMin() + ( m_dataField -> getXMax() - m_dataField -> getXMin() ) * x;
//    y = m_dataField -> getYMin() + ( m_dataField -> getYMax() - m_dataField -> getYMin() ) * y;
//    z = m_dataField -> getZMin() + ( m_dataField -> getZMax() - m_dataField -> getZMin() ) * z;
//    
//    sources[i].setCoordinates(x, y, z);
//  }

	m_dataField->getBound(m_xmin, m_xmax, m_ymin, m_ymax, m_zmin, m_zmax);
	m_Dim = (m_zmax - m_zmin < 1e-6) ? 2 : 3;
	Point *sources = NULL;
	if(2 == m_Dim)
	{
		//double x, y;
		//sources = new Point[m_sourceNumber];
	  //for(unsigned int i = 0; i < m_sourceNumber; i++)
	  //{
	  //  x = rand()/double(RAND_MAX);
	  //  y = rand()/double(RAND_MAX);
	  //  
	  //  x = m_dataField -> getXMin() + ( m_dataField -> getXMax() - m_dataField -> getXMin() ) * x;
	  //  y = m_dataField -> getYMin() + ( m_dataField -> getYMax() - m_dataField -> getYMin() ) * y;
	  //  
	  //  sources[i].setCoordinates(x, y, 0);
		//}

		size_t sz = (size_t)sqrt((double)m_sourceNumber) + 1;
		m_sourceNumber = sz * sz;
		sources = new Point[m_sourceNumber];
		double xStp = (m_xmax - m_xmin)/(sz+1);
		double yStp = (m_ymax - m_ymin)/(sz+1);
		size_t count = 0;
		for(size_t r = 1; r <= sz; r++)
		{
			for(size_t c = 1; c <= sz; c++)
			{
				sources[count++].setCoordinates( m_xmin + r * xStp,
																				 m_ymin + c * yStp,
																				 0);
			}
		}
	}
	else 
	{
		size_t sz = (size_t)pow(m_sourceNumber, 1.0/3) + 1;
		m_sourceNumber = sz * sz * sz;
		sources = new Point[m_sourceNumber];
		double xStp = (m_xmax - m_xmin)/(sz+1);
		double yStp = (m_ymax - m_ymin)/(sz+1);
		double zStp	= (m_zmax - m_zmin)/(sz+1);	
		size_t count = 0;
		for(size_t r = 1; r <= sz; r++)
		{
			for(size_t c = 1; c <= sz; c++)
			{
				for(size_t h = 1; h <= sz; h++)
				{
					sources[count++].setCoordinates( m_xmin + r * xStp,
																				   m_ymin + c * yStp,
																					 m_zmin + h * zStp
																					);
				}
			}
		}
	} 
	cout << "Sources number 	:	"	 << m_sourceNumber << endl;
	/******************* Set parameters and integrate ****************/
  visStreamline streamline(m_dataField);
	streamline.parseParameter(m_pfile);
  streamline.setIntegrationDouble();
  streamline.setSources(sources, m_sourceNumber);
  streamline.setParameters(m_stepLength, m_maxTime);

	time_t begin, end;
	begin = clock();
  m_streamlines = streamline.integrate();
	end = clock();
	cout << "Time used for streamline generation : " << (double)(end - begin)/CLOCKS_PER_SEC<< endl;

  delete []sources;
}

void visStreamlineClustering::computeStreamlinesInDomain()
{
	m_dataField->getBound(m_xmin, m_xmax, m_ymin, m_ymax, m_zmin, m_zmax);
	m_Dim = (m_zmax - m_zmin < 1e-6) ? 2 : 3;

	m_sourceNumber = m_Sources.size();
  Point *sources = new Point[m_sourceNumber];
	for(size_t i = 0; i < m_sourceNumber; i++)
	{ sources[i].setCoordinates(m_Sources[i].getX(), m_Sources[i].getY(), m_Sources[i].getZ()); }

	visStreamline streamline(m_dataField);
	streamline.parseParameter(m_pfile);
  streamline.setIntegrationDouble();
  streamline.setSources(sources, m_sourceNumber);
  streamline.setParameters(m_stepLength, m_maxTime);

	time_t begin, end;
	begin = clock();
  m_streamlines = streamline.integrate();
	end = clock();
	cout << "Time used for streamline generation : " << (double)(end - begin)/CLOCKS_PER_SEC<< endl;

	delete []sources;
}

void visStreamlineClustering::computeLength()
{
	/**
	 * Compute the streamline lenght infomation which will be stored int m_streamlineLenght.
	 */
	m_streamlineLength = new vector<double>[getStreamlinesNumber()];
	for(unsigned int i = 0; i < getStreamlinesNumber(); i++)
	{
		m_streamlineLength[i].push_back(0);
		for(unsigned int j = 1; j < getStreamlineSize(i); j++)
		{
			Point curr = getPointOnStreamline(i, j);
			Point prev = getPointOnStreamline(i, j-1);
			double dis = curr.Distance(prev);
			dis += m_streamlineLength[i][j-1];
			m_streamlineLength[i].push_back(dis);
		}
	}

	// Newly added..
	determineIntervalAndSeg();
}

void visStreamlineClustering::determineIntervalAndSeg()
{
	size_t dataPoint = m_dataField -> getPntNumberOfModel();

	if(2 == m_Dim)
	{
		m_Interval = (m_xmax-m_xmin)*(m_ymax-m_ymin)/dataPoint;
		m_Interval = sqrt(m_Interval);
	}
	else
	{
		m_Interval = (m_xmax-m_xmin)*(m_ymax-m_ymin)*(m_zmax-m_zmin)/dataPoint;
		m_Interval = pow( m_Interval, 1.0/3.0 );
	}
	double maxLen = -DBL_MAX;
	for(size_t i = 0; i < getStreamlinesNumber(); i++)
	{
		size_t sz = m_streamlineLength[i].size();
		double lineLen = m_streamlineLength[i][sz - 1];
		maxLen = maxLen > lineLen ? maxLen : lineLen;
	}
	m_seg = size_t((maxLen/m_Interval + 1));
//	cout << "xmax : xmin : " << m_xmax << ", " << m_xmin << endl;
//	cout << "ymax : ymin : " << m_ymax << ", " << m_ymin << endl;
//	cout << "zmax : zmin : " << m_zmax << ", " << m_zmin << endl;
//	cout << "data point number : " << dataPoint << endl;
	cout << "Max streamline len	:	" << maxLen << endl;
	cout << "streamline segment interval : " << m_Interval << endl;
//	cout << "streamline max len : " << maxLen << endl;
	cout << "streamline segments : " << m_seg << endl;
	cout << "SUGGESTED CLUSTERING THRESHOLD : " << m_Interval << endl;
}

Point visStreamlineClustering::getPointOnStreamline(unsigned int line, unsigned int pnt) const
{
  // Note that the first point is at the end of the streamline of the backward direction, while
  // the last point at the end of the streamline of the forward direction.

  // In the streamline of backward direction.
  unsigned int backward = line + getStreamlinesNumber();

  if( (*m_streamlines)[backward].size() > pnt )
  {
    return (*m_streamlines)[backward][ (*m_streamlines)[backward].size() - pnt - 1 ];
  }
  else
  {
  	pnt -= (*m_streamlines)[backward].size();
  	return (*m_streamlines)[line][pnt];
  }
  return 0;
}

void visStreamlineClustering::getPointOnStreamline
															( unsigned int line, 
															  const unsigned int segments, 
																unsigned int pnt,
																Point &p
															) const
{
//  /**
//	 *   lowIndex       pnt      highIndex
//	 *    *______________*_________*
//	 *    |<----ratio--->|
//	 *
//	 * The returned coordinates should be : lowIndex * (1-ratio) + highIndex * ratio
//	 */
//	size_t orignalSegments = getStreamlineSize(line) - 1;
//
//	/**
//	 * @Modification 2:
//	 * The original is:
//	 * 		double ratio = pnt / (double)segments * orignalSegments;
//	 * This really silly!
//	 */
//	double ratio = pnt / (double)segments * orignalSegments;
//	unsigned int lowIndex = (unsigned int)ratio;
//	unsigned int highIndex = lowIndex + 1;
//	ratio -= lowIndex;
//
//	Point *lowPnt = getPointOnStreamline(line, lowIndex);
//	Point *highPnt = getPointOnStreamline(line, highIndex);
//  p = highPnt->linearInterpolation(*lowPnt, ratio);
//
//
	if(0 == pnt)
	{
		p = getPointOnStreamline(line, 0);
		return;
	}
	if(segments == pnt)
	{
		p = getPointOnStreamline(line, getStreamlineSize(line)-1); 
		return; 
	}
	unsigned int sz = m_streamlineLength[line].size();
	double originalLen = m_streamlineLength[line][sz-1];
	double ratio = originalLen/segments*pnt;
	unsigned int highIndex = binarySearch(m_streamlineLength[line], ratio);
	unsigned int lowIndex = highIndex - 1;
	ratio -= m_streamlineLength[line][lowIndex];

	Point lowPnt = getPointOnStreamline(line, lowIndex);
	Point highPnt = getPointOnStreamline(line, highIndex);
  p = linearInterpolation(highPnt, lowPnt, ratio);
}

void visStreamlineClustering::getPointOnStreamline
															( const strmLine &line, 
																const vector<double> &lineLen,
															  const unsigned int segments, 
																unsigned int pnt,
																Point &p
															) const
{
	if(0 == pnt)
	{
		p = line[0];
		return;
	}
	if(pnt == segments)
	{
		p = line[line.size()-1]; 
		return; 
	}
	unsigned int sz = lineLen.size();
	double originalLen = lineLen[sz-1];
	double ratio = originalLen/segments*pnt;
	unsigned int highIndex = binarySearch(lineLen, ratio);
	unsigned int lowIndex = highIndex - 1;
	ratio -= lineLen[lowIndex];
	Point highPnt = line[highIndex];
	Point lowPnt = line[lowIndex];
  p = linearInterpolation(highPnt,  lowPnt, ratio); 

//  /**
//	 *   lowIndex       pnt      highIndex
//	 *    *______________*_________*
//	 *    |<----ratio--->|
//	 *
//	 * The returned coordinates should be : lowIndex * (1-ratio) + highIndex * ratio
//	 */
//	/**
//	 * @Modification 1
//	 * size_t orignalSegments = getStreamlineSize(line) - 1;
//	 */
//	size_t orignalSegments = line.size() - 1;
//
//	/**
//	 * @Modification 2:
//	 * The original is:
//	 * 		double ratio = pnt / (double)segments * orignalSegments;
//	 * This really silly!
//	 */
//	double ratio = pnt / (double)segments * orignalSegments;
//	unsigned int lowIndex = (unsigned int)ratio;
//	unsigned int highIndex = lowIndex + 1;
//	ratio -= lowIndex;
//
//	Point lowPnt = line[lowIndex];
//	Point highPnt = line[highIndex];
//  p = highPnt..linearInterpolation(lowPnt, ratio);
}

double visStreamlineClustering::EuDistance(const strmLine &line1, const strmLine &line2)
{
	double sim = 0;

	vector<double> len1, len2;
	computeLength(line1, len1);
	computeLength(line2, len2);
	double length1 = len1[len1.size() - 1], length2 = len2[len2.size() - 1];
	if( 0 == length1 || 0 == length2 ){	return DBL_MAX;}

	double *error = new double[m_seg + 1];
	for(unsigned int i = 0; i <= m_seg; i++)
	{
		Point pnt1, pnt2;
		getPointOnStreamline(line1, len1, m_seg, i, pnt1);
		getPointOnStreamline(line2, len2, m_seg, i, pnt2);
		error[i] = pnt1.Distance(pnt2);
	}
	double avg = 0, dist = 0;
	for(size_t i = 0; i <= m_seg; i++) {	avg += error[i]; }
	avg /= (1 + m_seg);
	for(size_t i = 0; i <= m_seg; i++) { dist += abs( avg - error[i] ); }
	dist /= (1 + m_seg);

	sim = /*0.65 * dist + 0.35 * */avg;
	delete []error;
	return sim;
}

double visStreamlineClustering::EuDistance(unsigned int line1, const strmLine &line2)
{
	double sim = 0;
	vector<double> len;
	computeLength(line2, len);
	if( 0 == len[len.size() - 1] || 0 == m_streamlineLength[line1][m_streamlineLength[line1].size() - 1])
	{ return DBL_MAX;	}

	double *error = new double[m_seg + 1];
	for(unsigned int i = 0; i <= m_seg; i++)
	{
		Point pnt1, pnt2;
		getPointOnStreamline(line1, m_seg, i, pnt1);
		getPointOnStreamline(line2, len, m_seg, i, pnt2);
		error[i] = pnt1.Distance(pnt2);
	}
	double avg = 0, dist = 0;
	for(size_t i = 0; i <= m_seg; i++){	avg += error[i]; }
	avg /= (1 + m_seg);
	for(size_t i = 0; i <= m_seg; i++){ dist += abs( avg - error[i] ); }
	dist /= (1 + m_seg);

	sim = /*0.6 * dist + 0.4 * */avg;
	delete []error;
	return sim;
}

double visStreamlineClustering::EuDistance(unsigned int line1, unsigned int line2)
{
	double sim = 0;
	Point pnt1, pnt2;
	if( 0 == m_streamlineLength[line1][m_streamlineLength[line1].size() - 1] || 
			0 == m_streamlineLength[line2][m_streamlineLength[line2].size() - 1])
	{ return DBL_MAX;	}

	double *error = new double[m_seg + 1];
	for(unsigned int i = 0; i <= m_seg; i++)
	{
		getPointOnStreamline(line1, m_seg, i, pnt1);
		getPointOnStreamline(line2, m_seg, i, pnt2);
		error[i] = pnt1.Distance(pnt2);
	}
	double avg = 0, dist = 0;
	for(size_t i = 0; i <= m_seg; i++){	avg += error[i]; }
	avg /= (1 + m_seg);
	for(size_t i = 0; i <= m_seg; i++){ dist += abs( avg - error[i] ); }
	dist /= (1 + m_seg);

	sim = /*0.6 * dist + 0.4 * */avg;
	delete []error;
	return sim;
}

double visStreamlineClustering::minDistance(size_t line1, size_t line2)
{
	double sim = 0;
	Point pnt1, pnt2;
	if( 0 == m_streamlineLength[line1][m_streamlineLength[line1].size() - 1] || 
			0 == m_streamlineLength[line2][m_streamlineLength[line2].size() - 1])
	{ return DBL_MAX;	}

	kdTree lineTree(m_seg + 1);
	for(size_t i = 0; i <= m_seg; i++)
	{
		Point pnt;
		getPointOnStreamline(line2, m_seg, i, pnt);
		lineTree.store( pnt.getX(), pnt.getY(), pnt.getZ() );
	}
	lineTree.treeBalance();

	const double disError = pow(m_xmin-m_xmax,2)+pow(m_ymin-m_ymax,2)+pow(m_zmin-m_zmax,2);
	for(unsigned int i = 0; i <= m_seg; i++)
	{
		Point pnt;
		getPointOnStreamline(line1, m_seg, i, pnt);
		nNearestNodes nNN(1);
		nNN.setDistError(disError);
		nNN.setSearchPnt( pnt.getX(), pnt.getY(), pnt.getZ() );
		lineTree.locateNodes(&nNN, 1);
		Point pntOnTree(nNN.index[1]->pos[0],nNN.index[1]->pos[1],nNN.index[1]->pos[2]);
		sim += pnt.Distance(pntOnTree);
	}
	sim /= (1 + m_seg);

	return sim;
}

double visStreamlineClustering::minDistance(strmLine &line1, strmLine &line2)
{
	double sim = 0;
	vector<double> len1, len2;
	computeLength(line1, len1);
	computeLength(line2, len2);
	double length1 = len1[len1.size() - 1], length2 = len2[len2.size() - 1];
	if( 0 == length1 || 0 == length2 ){	return DBL_MAX;}

	kdTree lineTree(m_seg + 1);
	for(size_t i = 0; i <= m_seg; i++)
	{
		Point pnt;
		getPointOnStreamline(line2, len2, m_seg, i, pnt);
		lineTree.store( pnt.getX(), pnt.getY(), pnt.getZ() );
	}
	lineTree.treeBalance();

	const double disError = pow(m_xmin-m_xmax,2)+pow(m_ymin-m_ymax,2)+pow(m_zmin-m_zmax,2);
	for(unsigned int i = 0; i <= m_seg; i++)
	{
		Point pnt;
		getPointOnStreamline(line1, len1, m_seg, i, pnt);
		nNearestNodes nNN(1);
		nNN.setDistError(disError);
		nNN.setSearchPnt( pnt.getX(), pnt.getY(), pnt.getZ() );
		lineTree.locateNodes(&nNN, 1);
		Point pntOnTree(nNN.index[1]->pos[0],nNN.index[1]->pos[1],nNN.index[1]->pos[2]);
		sim += pnt.Distance(pntOnTree);
	}
	sim /= (1 + m_seg);
	return sim;
}

double visStreamlineClustering::minDistance(unsigned int line1, strmLine &line2)
{
	double sim;
	cout << "Please chech ... wrong at !!" << endl;
	exit(-1);
	return sim;
}

double visStreamlineClustering::formalMetric(strmLine &line1, strmLine &line2)
{
	double sim = 0;
	return sim;
}

double visStreamlineClustering::formalMetric(unsigned int line1, strmLine &line2)
{
	double sim = 0;
	return sim;
}

//void visStreamlineClustering::computeDeSimilarity(pf1 pDeSimFun)
//{
//	m_relationMatrix = new complex_[getStreamlinesNumber() * getStreamlinesNumber()];
//	// Compute the desimilarity value
//	for(size_t row = 0; row < getStreamlinesNumber(); row++)
//	{
//		setSquareMat(m_relationMatrix, m_sourceNumber, row, row, 0);
//		for(size_t col = 0; col < row; col++)
//		{
//			double tmp = (this->*pDeSimFun)(row, col);
//			setSquareMat(m_relationMatrix, m_sourceNumber, row, col, tmp);
//			setSquareMat(m_relationMatrix, m_sourceNumber, col, row, tmp);
//		}
//	}
//}

double visStreamlineClustering::computeDeSimilarity(strmLine &line1, strmLine &line2, pf1 pDeSimFun)
{
	return (this->*pDeSimFun)(line1, line2);
}

double visStreamlineClustering::computeDeSimilarity(unsigned int line1, strmLine &line2, pf2 pDeSimFun)
{
	return (this->*pDeSimFun)(line1, line2);
}

void visStreamlineClustering::transformDeSim2SimMatrix()
{

}

/***********************************************Clustering***********************************/
double visStreamlineClustering::totalError()
{
	return 0;
}

void visStreamlineClustering::initClusters(unsigned int k, unsigned int iteration)
{
	if(k <= 0 || k > getStreamlinesNumber())
	{
		cout << "Invalide cluster nuber used at line " << __LINE__ << " of file " __FILE__ << endl;
		exit(-1);
	}
	m_itrt = iteration;
	m_NumClusters = k;
	/** 
	 * @Modification 4 : the m_ClusterCenters and m_ClassHaveObjVec are allocated memory of 
	 * m_NumClusters*(m_itrt+1) and m_NumClusters*(m_itrt+1) respectively to store the clustering 
	 * process.
	 * Note that the original choozen centers will be stored at the first m_NumClusters of the 
	 * m_ClusterCenters.
	 */
	m_ClusterCenters = new vector<Point>[m_NumClusters*(m_itrt+1)];
	m_ObjIsClassVec = new unsigned int[getStreamlinesNumber()];
	m_ClassHaveObjVec = new vector<unsigned int>[m_NumClusters*m_itrt];
	for(unsigned int i = 0; i < getStreamlinesNumber(); i++){ m_ObjIsClassVec[i] = -1;}
	if( !m_ClusterCenters || !m_ObjIsClassVec )
	{
		cout << "Memory out at line " << __LINE__ << " of file " << __FILE__ << endl;
		exit(-1);
	}

	unsigned int nonEmptyIndex = 0;
	for(unsigned int i = 0; i < m_NumClusters; i++)
	{
		// Initialize the initial cluster center as the first k NON-EMPTY ONES of the streamlines.
		while(nonEmptyIndex < getStreamlinesNumber() && 5 > getStreamlineSize(nonEmptyIndex))
		{	nonEmptyIndex++;}
		if(nonEmptyIndex >= getStreamlinesNumber())
		{
			cout << "Cluster iniitalization fails! Retry please..." << endl;
			exit(-1);
		}
		for(unsigned int j = 0; j < getStreamlineSize(nonEmptyIndex); j++)
		{
			Point pnt = getPointOnStreamline(nonEmptyIndex, j);
			setClusterCenter(0, i, pnt);
		}
		// Indicate that the i-th streamline belongs to class i.
		m_ObjIsClassVec[nonEmptyIndex] = i;
		nonEmptyIndex++;
	}
	//totalError();
}

void visStreamlineClustering::reshapeCluster2One(_SCHEME schm, unsigned int itrt)
{
	switch(schm)
	{
		case SIMPLEST:
		{
			for(unsigned int cls = 0; cls < m_NumClusters; cls++)
			{
				for(unsigned int numPnt = 0; numPnt < m_seg; numPnt++)
				{
					Point pnt, tmp;
					for(unsigned int obj = 0; obj < getClusterMembers(itrt, cls).size(); obj++)
					{
						getPointOnStreamline(getClusterMember(itrt, cls, obj), m_seg, numPnt, tmp);
						pnt += tmp;
					}
				  pnt /= getClusterMembers(itrt, cls).size();
					getClusterCenters(itrt, cls).push_back(pnt);
				}
			}	
			break;
		}// SIMPLEST
		case DIS_WEIGHTED:
		{
//			double ratio = 0.3;
			// If it's not zero-th iteration, then the following will access the previous result but not
			// the original(0-th). 
			for(unsigned int cls = 0; cls < m_NumClusters; cls++)
			{
				/**
				 * It's bad news that the cluster center can be destroied step by step!!
				 *
				 * Is it neccessary to re-allocate memory for the empty centers???
				 */
				if(0 == getClusterMembers(itrt, cls).size())continue;
				for(unsigned int numPnt = 0; numPnt < m_seg; numPnt++)
				{
					Point pnt, tmp;
					for(unsigned int obj = 0; obj < getClusterMembers(itrt, cls).size(); obj++)
					{
						// The index of line
						getPointOnStreamline(getClusterMember(itrt, cls, obj), m_seg, numPnt, tmp);
						pnt += tmp;
					}
				  pnt /= getClusterMembers(itrt, cls).size();
					setClusterCenter(itrt+1, cls, pnt);
				}
			}
			break;
		}//DIS_WEIGHTED
	default:
	{
		break;
	}
	}
}

unsigned int visStreamlineClustering::binarySearch(const vector<double> &data, double key)const
{
	unsigned int low = 0, high = data.size() - 1;
  unsigned int h = high;
  
  while(low <= high)
  {
    int mid = (high + low)/2;
    if( key == data[mid])
      return mid;
    
    if( data[mid] < key )
      low = mid + 1;
    else
      high = mid - 1;
  }

  // Search fails! Then return the position to insert the key.
  if( low <= h )
  {
    return low;
  }
  else
  {
    return h + 1;
  }

}

void visStreamlineClustering::computeLength(const strmLine &line, vector<double> &tmpLen)
{
	/**
	 * Compute the streamline lenght infomation which will be stored int m_tmpLen.
	 */
	tmpLen.clear();
	tmpLen.push_back(0);
	for(unsigned int j = 1; j < line.size(); j++)
	{
		Point prest = line[j];
		Point pre = line[j-1];
		double dis = prest.Distance(pre);
		dis += tmpLen[j-1];
		tmpLen.push_back(dis);
	}
}

}// Namespace visualization
